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I. INTRODUCTION 

In the preceding paper |jj , the problem of finding ground state energies and configurations for a Frenkel-Kontorova 
model in a periodic potential formed by parabolic segments of identical (positive) curvature, was reduced to that of 
minimizing a certain convex function over a finite simplex. While various aspects of the corresponding phase diagram 
in the N = 2 case can be worked out in a relatively straightforward manner ^-Q] , minimizing the convex function for 
larger values of N represents a non-trivial problem in numerical analysis. The basic reason is that the convex function 
of interest does not have continuous derivatives, and in the case of an irrational winding number, it possesses a dense 
. set of singularities. Hence standard gradient methods run into difficulties. 

The approach we employ is based upon the concepts of subdifferential and subgradient from the theory of convex 
functions ||, as explained in Sec. ||. The algorithm itself, described in detail in Sec. [II, is motivated by a physical 



■ model involving quasiparticles, Sec. pi A. An essential part of the procedure is a standard linear programming 



procedure which we have simplified and adapted to the problem at hand, App. A, so as to speed it up substantially. 
As a result, ground states for N on the order of 100 arc readily calculated, and larger values of N are accessible with, 



of course, a longer running time; see Sec. Ill D 



II. MINIMIZATION USING SUBDIFFERENTIALS 



As in part I, with some minor changes in notation, we assume the energy per particle can be written in the form 
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where 



with ho equal to zero, 



and 



e = e + ei, (1) 



e = 5^At i At fc a(C«). ( 2 ) 

j<k 

N-l 

e 1 =h-ip = -f)-C = -^2fji&, (3) 

i=0 



Cfcj = -o* = Cfc - o = ( 4 ) 

i=2+l 



g(V,) = g(-^) = 0(i-v)= £ W^ 1 ^)- (5) 
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The boldface letters denote component vectors, for example, 

C = (&,&,». Gv-i). (6) 

A constant term independent of C has been omitted from (^), and bars have been added to h and rj in (||) to distinguish 
them from quantities which we shall define later. Note that 

N-l 

$> = ° (7) 

i=0 

because ^ = /ij+i — /ij, and /ij is periodic in j, with period AT. 

Our task is to find the £, or equivalently tp, which minimize e for a given i), or equivalently h. The function G(ip) 
is convex on the interval < ij> < 1. If its derivative 

= d^/dV = -B(-^) = B(l + i>) (8) 
were a continuous function, the minimum would satisfy the equation: 

V = V (9) 
obtained by differentiating (Q) , where the components of T] are given by 

Vk = Vkj, (10) 

r/ kj = AtfeAtj (3 kj = -7]jk, (11) 

kj = B(Gy) = -p jk . (12) 

But in fact B(ip) has lots of discontinuities, see Fig. [y, and therefore we need some way to interpret the (formal) 
solution (^|) in this case. For this purpose it is convenient to employ the concepts of subgradient and subdifferential 
as defined by Rockafellar ||. Suppose that F(Q is a real- valued function (which need not be convex) defined on some 
domain in M. N . We shall say that r/ is a subgradient of F at the point C = C provided 

F(C)>F(C)+v(<:-0 (13) 
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for all £ where F is defined. The collection of all rj values for which this inequality is satisfied for a given £ is easily 
shown to be a convex subset of R*, and is called the subdifferential of F at C, denoted by dF{C,). 

Given this definition, it is easy to show that e in ([!]) has a minimum at £ = £ if and only if rj is an element of 
9eo(C)j the subdifferential of eo at £. In addition, the subdifferential of a sum (F\ + F% + • - •) of convex functions is 
the sum of the subdifferentials of the individual functions Q , understood as sums of sets of R w ; thus: 

A + £? = {c + C':Ce AC'eB}, (14) 

together with its obvious generalization to the sum of three or more sets. 



These observations provide the key for interpreting equations (10) to (|l2|). At some point ip' where B(ip) is 



discontinuous, the subdifferential of Q{ip), for tp in the range < il> < 1, consists of all points (3 lying in the interval 

B-W) <P<B + W) (15) 

where B-(ip') and B+(ip') are the left and right derivatives of G(il>) at ip' , that is, the bottom and top of the 
discontinuity in the graph of B. More generally, for k > j we interpret f3kj in ( |l2| ) as any point in the interval 

MGy) < fly < B+(Cy)> (16) 

where the lower limit is set equal to — oo if Qkj = 0, and the upper limit is +oo if £jy = 1, as a consequence of the 
constraints 

Co <Cl <C2 < ■■■ <0v-i <Co + l- (17) 

For j < fc, we define f3jk = —@kj- Of course, if B{ip) is continuous at ip — Cfch then fluj is the single point B(Cki)- 

Consequently, 9eo(C) is simply the collection of all r/ obtained using (jlCf), where for each k > j, (3kj in ( M ) is 
allowed to vary over the interval ([H|), and for j < fc, f3jk = —fikj- Notice that this means that rjjk + r\kj is zero, and 
therefore 



N-l 

E 1 



= (18) 



for any r/ in the subdifferential deo(C), which corresponds to (0). Note that the f3kj are allowed to vary independently, 
aside from the restriction f3jk — —fikj- As there are thus N(N — l)/2 independent variables, some of which may be 
constant because they do not correspond to discontinuities of B, <9eo(C) is, in general, a fairly complicated polyhedron, 
of dimension less than or equal to N — 1, the dimension of the space in M. N satisfying the constraint jl^). 

In the case N = 3, the subdifferentials deo(C) are closed sets, either hexagons, lines, or points, depending upon the 
value of C- Some of the lines and hexagons extend to infinity. A hexagon occurs provided 

(lO^fllW-Vi, C20=M2W-^ 2 , (19) 

where fj,±, /i2, v\, and i>2 are integers, in which case £io and £20, as well as £21 = C20 ~ CiOj are at discontinuities of B. 
A line occurs when B is discontinuous at one of the three values ( w , ( 2 q, or £ 2 i) but not at the other two, and <9e (C) 
is a point if B is continuous at all three values. If w is irrational, the discontinuities of B(ip) are a dense set in tp, and 
consequently the subdifferentials of eo for different C have no points in common. If w is rational, adjacent hexagons 
overlap at their common edges and vertices, and each edge and each vertex is itself a subdifferential of eo for a range 
of £ values. For both rational and irrational w, the hexagons cover the entire plane satisfying the constraint (|lq), 
with the exception, when w is irrational, of a set of zero measure. A similar comment applies to larger values of N, 
and therefore in numerical studies it suffices to consider the N — 1 dimensional polyhedra obtained when every £kj 
falls on a discontinuity of B. 

It is sometimes helpful to think of the collection of subdifferentials 9eo(C) as C varies as generated by placing a set 
of N(N — l)/2 points on the graph of B(ip) at positions (Ck±, Pkj)- Note that the £fcj cannot be varied independently, 
as they are determined by a set of N — 1 parameters, see ([|). However, each of the iV(jV — l)/2 points on the graph 
can be moved independently in the vertical direction, as long as it is on a discontinuity of B, to form the collection of 
(3^ j values which generate the subdifferential for a fixed £. 



3 



III. NUMERICAL PROCEDURE 



A. Introduction 

The problem of finding the £ which minimizes e, (|l]) , for a given 77 is equivalent, as noted in Sec. || above, to 
finding the £ such that f) falls in the subdifferential <9eo(C). Furthermore, the N — 1 dimensional polyhedra which 
arise when all the (kj fall at the discontinuities of B(ip) fill up the relevant N — 1 dimensional hyperplane (jlj) except 
for a set of measure zero, and as we assume that r\ is only specified with some limited numerical precision, we can in 
practice limit ourselves to a consideration of such polyhedra. 

The general idea of the algorithm is as follows. Starting from some C with all (kj at discontinuities of B, test 
whether the target f) lies inside deo((). If it does, the problem has been solved. If it does not, use the information 
obtained from the test in order to choose a new ( closer to the desired value, and repeat the test. The test itself, 
steps 3 and 5 in the algorithm as summarized below, involves a linear optimization procedure with an execution time 
which (typically) varies as N 3 , which is relatively expensive when N is large. Consequently, the test is preceded in 
our algorithm by various steps whose aim is to provide, with a relatively small number of operations, a value of ( 
close to the final solution. 

To begin with, we replace the actual B(ip) with an approximate, piecewise constant function B*{tp) which has 
discontinuities at the points 

ip = (iw-v, (20) 

where fi and v are integers, and 

\lA < M (21) 

for some finite bound M, which can be increased later if necessary. This is a sensible procedure, because the size of 
the discontinuities decreases exponentially with Between two successive discontinuities ip 1 and i/>", the function 
B* is defined to be a constant lying halfway between B+(ip') and see Fig. |l|. Consequently, the discontinuities 

of B* are somewhat larger than those of the exact B, and ( |l6|) is replaced by: 

B*_(( kj )<0 kj <B* + (( kj ). (22) 

Note that if 

w = p/q (23) 

is a rational number, with p and q relatively prime positive integers, the discontinuities of B, (pOj), are the points 

if, = s/q, (24) 

where s is any integer (and may have factors in common with q) . The definition of B* is the same as before; though 
it should be noted that the discontinuity interval of B at a point ( p4| ) is made up of contributions from an infinite 
number of discontinuities from derivatives of terms on the right side of (||). Of course, if M in ( pil| ) is equal to q — 1 
(or larger), B* and B are identical, and step 5 can be eliminated from the algorithm described below. 

In order to motivate the initial steps in the algorithm, it is helpful to think of Coj Ci> • • ■ as representing the positions 
of a set of N quasiparticles located on a circle of unit circumference, Fig. ^, and subjected to two kinds of forces, 
corresponding to eo and e\, thought of as potential energies. In this picture, fjk represents an external force exerted 
on particle k, and rjkj, (0), the force which particle k exerts on particle j. The minimization condition (||) can then 
be interpreted as stating that, for every fc, the external force exerted on particle k is equal to the sum rjk of the forces 
which it exerts on the other particles, which is the same as saying that the net force on particle fc is zero. The pair 
force r/kj is constant as long as (kj is not at a discontinuity of B* , while if it is at such a discontinuity, it can take any 
value, see (|ll|), corresponding to (3kj in the range (^). In addition, there is a hard core interaction which prevents 
two quasiparticles from passing through each other, and ensures that the inequalities (|l7| ) are satisfied. Note that 
there can very well be solutions to the minimization problem in which some of these inequalities are equalities. If, for 
example, (2 = (3, then ^32 can be very large and negative, see the remarks following (|l6|), corresponding to the fact 
that the hard core allows particle 3 to exert a very large (negative) force on particle 2. 

The initial steps of the algorithm consist of a number of horizontal and vertical shifts; the terminology comes from 
the picture of points on the graph of B, Sec. O. A horizontal shift is a change in the positions of the quasiparticles, and 
thus the (kj ' with the (3kj (and thus the r/kj ) held fixed, while a vertical shift is a change in the set of (3kj values with 
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the quasiparticle positions, and thus the Qkj , held fixed. In addition, we shall make use of the concept of a cluster, 
which means a collection of quasiparticles, with their labels belonging to an index set J containing |J| members, 
with the property that the collection is connected by a set of "pair bonds" (kj), k and j members of J, with £jy at 
one of the discontinuities of B*. For example, if (21 an d C42 fall on discontinuities of B*, then it is possible, but not 
necessary, to define a cluster J = {1,2,4} of \ J\ = 3 quasiparticles, as in Fig. B. We shall always think of the entire 
collection of quasiparticles as divided up among a set of mutually disjoint clusters, where a quasiparticle which does 
not belong to a larger cluster constitutes its own cluster containing only one element. If all the quasiparticles belong 
to a single cluster of | J| = N elements, we shall call this a complete cluster. In any horizontal shift, the clusters are 
moved rigidly, in the sense that Qk does not change if j and k belong to the same cluster. (This must obviously be 
the case whenever the cluster is linked together by bonds for which the ftkj fall in the interiors of the corresponding 
discontinuity intervals (|22).) Conversely, a vertical shift is always applied to a single cluster. 



B. Summary of the Algorithm 



The algorithm for finding a minimum consists of the following steps, details of which are given below in Sec. Ill C . 

0. Initialization: 

Choose an initial approximate B* by, for example, setting M = 4 in (pl|), and some initial values for the Q satisfying 
(|l7|), with a set of clusters specified (e.g., each quasiparticle might belong to its own cluster). 

1. Horizontal shift I: 

Calculate a "velocity" for each cluster, and use this to carry out a horizontal shift until for the first time some £fcj 
for j and k belonging to different clusters reaches a discontinuity of B* , in which case we shall say that these two 
clusters have "collided" to form a temporary combined cluster. Go to step 2. 

2. Vertical shift: 

Carry out a vertical shift on the temporary combined cluster following the prescription given in ( |27| ) below and 
in the remarks which follow, and apply the test which is described there. If the result of the test is negative, the 
combined cluster is rejected, the collection of clusters is defined to be the same as it was before the collision, and the 
algorithm returns to step 1. If the result of the test is positive, the temporary combined cluster becomes permanent, 
and is considered part of the collection of clusters for the next step in the algorithm. If this cluster is complete, go to 
step 3; if not, return to step 1. 

3. Linear optimization I: 



With all the quasiparticles in a single cluster, apply linear optimization, as discussed in Sec. Ill C below, to produce 
a vertical shift which maximizes a non-negative parameter A. If A > 1, then the target f) is inside the polyhedron 
<9eo(C) associated with the current £ in the approximation in which B has been replaced by B*; go to step 5. If A < 1, 
then 77 is not inside the polyhedron; go to step 4. 

4. Horizontal shift II: 

Use the (3k j resulting from the linear optimization in 3 to divide the collection of quasiparticles into two clusters, 
which undergo a horizontal shift relative to each other until they collide (as in step 1) to form a new, combined cluster 
which is a complete cluster. Return to step 3. 

5. Linear optimization II: 

Repeat the linear optimization step 4, but with each /3kj now restricted to the corresponding exact interval (|ltl). 
If, however, £jy falls at a point where B* , unlike B, has no discontinuity, then the corresponding (3kj is placed at 
the center of the corresponding discontinuity of £>, and is treated as a constant, not a variable, during the linear 
optimization. If the optimization yields A > 1, the current £ is the desired solution to the minimization problem, and 
the algorithm stops. If A is less than 1, B* is replaced by another approximation to B constructed in the same way, 
but using a larger value of M in (|2l]) . The current £ values are changed by very small amounts so that none of the 
C,kj fah at discontinuities of the new £>*, and the algorithm returns to step 1. 



C. Details of the Algorithm 

The explanations given below are numbered in the same way as the steps in the preceding summary. 
1. Given a set of rj^j values, the net force on the fc'th quasiparticle is 

fjk - Vk = fjk - Vki- ( 25 ) 
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Were the force given by a continuous function, it would be possible to find the energy minimum by assigning to 
each quasiparticle a velocity proportional to the force acting on it, and then solving the resulting dynamics. What 
is actually done in the algorithm is to assign to each cluster a velocity given by the total force acting on all the 
quasiparticles in the cluster divided by the number of particles in the cluster, which is the average force per particle: 

feeJ 

It is only the relative cluster velocities which are of interest in determining the horizontal shift; adding the same 
constant to every vj will make no difference, and one can arrange (for example) that the cluster containing £o remains 
fixed. The clusters are then shifted by amounts proportional to their respective velocities until the first "collision" 
occurs, in the sense that for k in one cluster and j in another reaches a discontinuity of 23*. 

2. Once the temporary, combined cluster J c has been formed, a vertical shift is applied to the rj k j for k and j in 
J c . This is done by first calculating a preliminary value for the shift of the rj k j or f3 k j values from the formula 

Srjkj = AtfcAtj 5(3 kj = [(fj k - rjk) - (Vj ~ Vj)]/\Jc\- ( 27 ) 

The motivation for this choice is the following. If all the quasiparticles were in a single cluster, \J C \ = N, the change 
( p7| ) would result in a new set of pair forces 

rf kj = At k At, ki = Vkj + 5r) kj (28) 

with the property that 

that is, one would have solved the minimization problem. With \J C \ < AT, the result would, instead, be to make the 
difference fjk — rj' k , the sum of the forces acting on quasiparticle k, independent of k for all k in J c , and to minimize 

£(%-^) 2 (30) 

fceJ c 

as much as is possible by changing only the pair interactions rj k j inside the cluster J c . 

However, formula fl^ ) does not take account of the possibility that the in ( p8| ) might lie outside the interval 
( f22| ) determined by 23* . When this is the case, the new /3' k - is placed at whichever end of the discontinuity interval lies 
closest to the value given by (|8|). The test for rejecting or retaining the combined cluster J c is then the following. If 
for each of the pairs k and j for which k belongs to one of the clusters involved in the collision and j to the other, the 
new /3L is at one of the ends of the interval (^2|), the combined cluster J c is rejected, whereas if at least one of these 
values falls in the interior of the corresponding interval, J c is accepted. Note that whether the cluster J c is accepted 
or rejected, the new values of and thus the corresponding rf k j, produced in the vertical shift are retained when 
going on to the next step of the algorithm, which is either step 1 or, in the case in which J c is accepted and \J C \ = N, 
step 3. 

The algorithm would still function correctly if a combined cluster were never rejected. However, this would mean 
having to apply linear optimization, step 3, more often, and would result in a slower computation. The process of 
allowing clusters to move relative to each other past discontinuities which represent relatively small changes compared 
to the large forces representing a situation "far from equilibrium" helps to achieve a better preliminary value of C 
before going on to step 3. 

3, 4. The linear optimization step is basically a test to see whether the target f) lies inside the polyhedron 
representing <9eo(C) m t ne approximation in which 23 is replaced by 23*. The idea is to begin at a particular vertex rf 
of the polyhedron, Fig. ||, and draw a straight line from rf to the target f). Points along this line are of the form 

V = r, c + \(r l -r l c ), (31) 

where A is a number between and 1. The linear optimization procedure, the details of which are given in the 
appendix, determines the largest value of A for which a point of the form ( fjl| ) lies inside or on the boundary of 
the polyhedron. If this value is less than 1, as in Fig. ||, the target lies outside the polyhedron, and the point ( |3l| ) 
determined by the maximum value of A specifies a facet of the polyhedron lying in the direction of the target, as viewed 
from the starting vertex rf. This facet is generated by letting (N — 2) of the f3 k j vary over their entire discontinuity 
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intervals, while the remaining f3kj are fixed either at the top or at the bottom of their discontinuity intervals. The 
Qkj corresponding to the former are "rigid" in the sense that they cannot be altered by a horizontal shift (which, by 
definition, must leave the flkj unchanged), and one can identify two clusters of quasiparticles, each one connected by 
such rigid bonds. 

Once these two clusters have been identified, they can be shifted relative to each other, in a direction which is 
obvious, until they collide at the first discontinuity of B*. This collision results in a new, complete cluster, and the 
corresponding <9eo(C) is a polyhedron adjacent to the one considered earlier, and shares with it the facet which was 
identified in the previous linear optimization step. 

Note that this new cluster is accepted without carrying out the test used in step 2 of the algorithm. Also, in the 
unlikely event that the maximum A corresponds to the intersection of two or more facets of the polyhedron, the actual 
optimization algorithm described in the appendix will, in effect, "choose" one of these facets, and thus the polyhedron 
adjacent to it in the direction of the target r). 

5. If the optimization carried out in step 3 yields A > 1, the target f) lies inside the polyhedron generated by the 
discontinuities of B* , but it may or may not lie inside the corresponding polyhedron generated by restricting the /3kj 
to lie in the exact interval ( |16|) for the corresponding discontinuity. To test whether this is the case, one repeats 
the linear optimization technique of step 3, but now starting with each (3kj at its exact maximum possible value, and 
constrained to be greater than or equal to its exact minimum possible value, with the exception of those flkj for which 
Qkj does not fall at a discontinuity of B* , which are assigned fixed values at the center of the appropriate (exact) 
discontinuity intervals of B. One could, of course, make all of the variable during the linear optimization process. 
However, as the discontinuity intervals in B which are not in B* are, by construction, relatively small, the main effect 
of using a larger set of variables would be to slow down the linear optimization without much hope of actually finding 
a solution with A > 1 when the restricted search yields one with A < 1. Note that if such a restriction does result in 
overlooking a A > 1 solution, this solution, or one equivalent to it, will nevertheless be found later when the number 
of discontinuities of B* is increased. 

If the maximum value of A obtained by using linear optimization with the exact discontinuity intervals of B is 
still greater than or equal to 1, then the current £, and thus the current set of represents an actual solution to 
the minimization problem, and this does not depend upon the approximations used in constructing B* , because the 
resulting f3kj values all fall within the range where B* is identical to B. 

If, on the other hand, one finds that A is less than 1, this means that the current £ is not a solution to the true 
minimization problem; instead, it is as good as one can do using the approximate B* . To do better, it is necessary to 
increase the number of discontinuities. Our procedure at this point is to throw away all the information associated 
with the Pkj values obtained in the immediately preceding step of linear optimization, and simply start over again at 
step 1 using the current C- One might be able to improve the algorithm in this respect, but since the initial steps of 
the algorithm are relatively fast, it does not seem likely that one would obtain a significant increase in speed. 

D. Implementation and Performance 

The program we constructed to implement the algorithm described above was tested in the following way. We chose 
a winding number equal to the inverse golden mean (0.618. . . ) and a value of k, see [[!), of 0.6, resulting in a set of 
60 pairs {tp and 1 — "0) of discontinuities of B in the interval < tp < 1 with a magnitude greater than a resolution 
of 10 -20 . (Some tests used alternative values for k, 0.1, 0.5, and 1.0, for which there are 146, 66, and 47 pairs of 
discontinuities, respectively, exceeding this resolution.) The parabolas were assumed to be equally spaced, with At; 
independent of I. Then we employed the following "inverse strategy". With N fixed, random values of Q, lying on 
the full set of discontinuities of B were chosen, subject to the constraints (O) , and values rjkj inside the discontinuity 
intervals were also chosen randomly, thus defining — see (|Io|), and (|3j) — f) for a model of N parabolas with the 
solution to its energy minimization problem already known. The algorithm was then applied to this model starting 
at the initialization step 0, with a (different) random collection of Ci, and a choice of M (the number of pairs of 
discontinuities in the approximate £>*) to search for the correct solution. 

Note that the running time increases linearly with M. But since as M increases, the size of the discontinuities 
in B* is decreasing, the probability that a random choice of f) will actually require a larger value of M goes to zero 
exponentially with increasing M. We found that using an initial value of M = 2 when N is small saves a lot of 
running time, and in almost all cases M = 10 was sufficient to find a ground state with N up to 100. 

To determine the N dependence of the running time, we generated and timed 10 distinct potentials for each N in 
an increasing series up to TV = 115. Using an HP 9000 model 735 workstation with 64 MB of RAM and a CPU with 
20 MFLOPs, the average time in seconds required to find a solution was approximately 10~ 5 7V 4 , or a quarter of an 
hour for N = 100. The time required for linear optimization varies as TV 3 , but as iV becomes larger, "quasiparticle 
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dynamics" takes up a larger fraction of the time: 60-70% for iV 
in particular cases for N as large as 200. 



100. The algorithm found the correct ground state 
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APPENDIX A: LINEAR OPTIMIZATION 



The linear optimization procedure used in steps 3 and 5 of the algorithm functions in the following way. Written 
out in terms of components, Eq.(|3~l|) has the form: 



Vk = v c k + Hm -v c k)=Yl At k^tj/3kj 

j^k 



(Al) 



for k = 1, 2, . . . , N — 1. The equation for k 
fikj for k > j can be written as 



can be omitted because of the constraints (|7|) and (Tt8f). In turn, the 



= B;(Cy) + Xkj[BL(Cki) - BT + (Ckj)], 
where the N(N ~ l)/2 variables Xkj take values on the closed interval 

< x kj < 1. 



(A2) 



(A3) 

In step 5 of the algorithm, B* in ( |A2| ) should be replaced with the exact function B. The vertex rf of the polyhedron 
is defined as the value of rj when A = 0, and is obtained by setting all the Xkj in (A3) equal to zero, and inserting the 
resulting /3kj in the right side of JA1| ). It is, thus, a known quantity, as is also the t arge t f). By contrast, the xuj and 
A are considered as variables subject to the N — 1 equalities obtained by inserting (A3) in (Al), and the inequalities 
(A3), along with A > 0. If we define the objective function to be equal to A, maximizing A subject to these constraints 
becomes a standard problem in the theory of linear optimization. For a compact but clear description, see 0. The 
usual approach is to replace the second set of inequalities in (A3) by the requirement that the "slack" variables 



Vkj 1 %kj 



(A4) 



be non-negative. However, in order to minimize the number of variables, we have modified the formulation in 0, as 
discussed below, so that at any point in the calculation, either ajjy or y^j, but not both, appears as a variable in the 
tableau. 

The first task is to construct a tableau in restricted normal form, in which the N 



1 equalities (Al), expressed in 
1 of the Xkj variables are expressed as 
equation is missing 



terms of the Xkj using (A2), are transformed into expressions in which N 
linear functions of the remaining x's and A. To see how this is done, note that since the fc 
from the set (Al), X ko o ccurs only in the fc'th equation, and in none of the others. Therefore, if, for each k > 0, the 
coefficient of Xko in (A2) is non-zero, we can rewrite the N — 1 equalities in the form 



XkO 



(A5) 



for k = 1, 2, . . . , N — 1, with no xjq appearing anywhere on th e rig ht side of any of these equations. 

It may, however, happen that the coefficient of some Xko in ( A2) vanishes, because B*(tf>) does not have a disconti- 
nuity at ip = Cfco- In this case the following strategy can be employed. Construct a graph in which the N quasiparticles 
are vertices, and the y are connected together in one cluster using N — 1 edges (k,j) chosen so that for each edge the 
coefficient of xj-j in (A2) does not vanish. Such a graph is a tree with no closed loops, because a connected graph 
of TV vertices which includes a loop will have at least N edges. O ne t hen chooses some vertex m which has only one 
edge connecting it to another vertex n, and solves the equation (|Al[ ) with k = m for x mn (or x nm ) in terms of the 
other it's. After this, eliminate vertex m and the edge (mn) from the graph, and repeat the process. At each stage 
one finds (it is useful to work out an example) that the only x's appearing on the right side of the equation either 
do not correspond to edges of the original graph, or else to edges which have already been eliminated. The latter, 
however, can be expressed (recursively) as functions of the x's corresponding to edges which were not present in the 
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original graph, thus leading to a set of N — 1 equalities of a form similar to ( A.5). Note that the Xkj whose coefficients 
are zero in (A2) can be completely eliminated in constructing the linear optimization tableau, as they play no role. 

Once the tableau has been constructed for the restricted normal form, the optimization is carried out as indicated 
in 0|, with the following difference. When considering increasing the variable associated with a particular column, it 
is necessary to take account of both positive and neg ativ e entries in the column in order to test whether one of the 
left will reach either the upper or lower limit in (|A3| ) If the choice is to set a left variable Xkj equal to 1, it is 
replaced as a right variable by the corresponding ykj , which is 0, and the tableau is appropriately transformed. 
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FIG. 1. The exact B(ip) (solid lines) for a winding number equal to the inverse golden mean (.618. . . ) and n — 0.1 is shown 
by solid lines, and the approximation B*(ip), which retains only the largest discontinuities (see text) by dashed lines. 



FIG. 2. An example of N = 6 quasiparticles lying on a circle of unit circumference and partially coupled to a cluster with 
index set J = {1, 2,4}. 



FIG. 3. Schematic representation of the linear optimization procedure, see (J3l|) . Points 77 with A less than 0.6 lie inside the 
polyhedron (in this case a polygon), so that the target f) lies outside it, and the open circle corresponding to the maximum 
value of A indicates the starting point for the next step of the algorithm. 
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